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Abstract. Computation of high frequency solutions to wave equations is important 
in many applications, and notoriously difficult in resolving wave oscillations. Gaussian 
beams are asymptotically valid high frequency solutions concentrated on a single curve 
through the physical domain, and superposition of Gaussian beams provides a powerful 
tool to generate more general high frequency solutions to PDEs. An alternative way to 
compute Gaussian beam components such as phase, amplitude and Hessian of the phase, 
is to capture them in phase space by solving Liouville type equations on uniform grids. In 
this work we review and extend recent constructions of asymptotic high frequency wave 
fields from computations in phase space. We give a new level set method of computing 
the Hessian and higher derivatives of the phase. Moreover, we prove that the k th order 

phase space based Gaussian beam superposition converges to the original wave field in 

2 h — n. 

L at the rate of £2 4 in dimension n. 
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1. Introduction 
In this paper we consider the following equation 

(l.i) PV = o, (x,t)eR n x 
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where P = —iedt + H(x, —ied x ) is a linear differential operator with a real principal symbol 
r + H(x,p), subject to the highly oscillatory initial data 

(1.2) i/>{x, 0) = := A m (x)e iS ^/*, 

where Ai n 6 C^°(R n ) and S{ n G C°°(IR n ). The canonical example is the semi-classical 
Schrodinger equation with the Hamiltonian H(x,p) = ^\p\ 2 + V ex t(x), where V ex t(x) is a 
given external potential. The small parameter e represents the fast space and time scale 
introduced in the equation, as well as the typical wave length of oscillations of the initial 
data. Propagation of oscillations of wave length 0(e) causes mathematical and numerical 
challenges in solving the problem. In this article we are interested in the construction 
of globally valid asymptotic wave fields and the analysis of their convergence to the true 
solutions of the initial value problem. 

Geometric optics, also known as the WKB method or ray-tracing, when applied to 
model high frequency wave propagation problems such as (II. ip leads to the WKB-type 
system for both phase and amplitude. The phase is governed by the Hamilton-Jacobi 
equation 

(1.3) d t S + H(x,V x S) = 0, xeR n , t > 0. 

Solving this equation using the method of characteristics can lead to singularities which 
invalidate the approximation. In general, this breakdown occurs when the density of rays 
becomes infinite. This corresponds to the formation of a caustic where geometric optics 
incorrectly predicts that the amplitude of the solution is infinite. The consideration of 
these difficulties, beginning with Keller [13j . and Maslov [23], led to the development of 
the theory of Fourier integral operators, e.g., as given by Hormander [7J. 

A closely related alternative to Fourier integral operators is the construction of approxi- 
mations based on Gaussian beams. Gaussian beams are asymptotic solutions concentrated 
on classical trajectories for the Hamiltonian H(x,p), and they have a history going back 
to at least the late 1960's. They were initially used to study resonances in lasers [I], and 
later to obtain results on the propagation of singularities in solutions of PDE's [8] and 
[28j . At present there is a considerable interest in using superpositions of beams to resolve 
high frequency waves near caustics. This goes back to the geophysical applications in [2] 
and [B]. Recent works in this direction include [32] on "gravity waves", |14] and [12j on 
the semi-classical Schrodinger equation, and [33] and [25] on acoustic wave equations. 

An alternative to the standard WKB method is to use multi- valued solutions, {Si(t, x)}fL 
to (jl.3p corresponding to crossing waves [31]. This is in sharp contrast to the notion of 
viscosity solution [4] adopted when (11.3P arises in other applications. In the last decade 
a considerable amount of work has been done to capture multi-valued phases associ- 
ated to the WKB system numerically; we refer to review articles [5j [30] and references 
therein. Recently the level set method has been developed to resolve multi-valuedness 
of involved quantities in phase space as well as to compute physical observables, e.g. 
[261 13 QH EH M 13 OH [El Q21 EB EI]. The key idea, for instance in [3 E7J , is to repre- 
sent characteristic trajectories by common zero sets of some implicit level set functions, 
and evolve all relevant quantities in phase space, see the review [18]. Another phase space 
based approach is the use of the Wigner transformation [M] to map solutions of the under- 
lying wave equation to functions on phase space. See [161 [231 E21 [31] for the application of 
the Wigner transformation to the semi-classical Schrodinger equation. These phase space 
based approaches are extremely useful since they unfold 'caustics'. However, at caustics, 
neither gives correct prediction for the amplitude. This has led to an effort to combine the 
accuracy of beam superpositions at caustics with the level set method of computation in 
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phase space. This approach with attention to the accuracy of the resulting approximations 
is the subject of this paper. 

In this article we have two objectives: 

(i) to present the construction of beam superpositions by level set methods ; 

(ii) to estimate the error between the exact wave field and the asymptotic ones. 

The construction for (i) is based on Gaussian beams in physical space, but it is carried 
out by solving inhomogeneous Liouville equations in phase space as in [15], [13] and [12] . 
The result is no longer a superposition of asymptotic solutions to the wave equation (1.1)! 
Since it can be written as a superposition of standard Gaussian beams composed with a 
time-dependent symplectic change of variables, the superposition over phase space is still 
an asymptotic solution, as was pointed out in [15j . Here we consider superpositions over 
subdomains moving with the Hamiltonian flow, and show directly that they are asymptotic 
solutions without reference to standard Gaussian beams. For (ii) we use the well-posedness 
theory for (1.1), i.e. the continuous dependence of solutions of Pip = / on their initial 
data and /. Thus, the sources of error in the Gaussian beam superposition for an initial 
value problem are the error in approximating the initial data and the error in solving the 
PDE. 

To be specific, our asymptotic solution is expressed as 

(1.4) ,p% y) = Z(n, e) / 4> PG B(t, y, X)dX, 

JQ(t) 

where X = (x,p) denotes variables in phase space M 2n , Q(0) is the domain where we 
initialize Gaussian beams from the given data, and Q(t) = X(t, 17(0)) is the image of fi(0) 
under the Hamiltonian flow. Here ipPGB(t,y, X) is the phase space based Gaussian beam 
Ansatz, and Z(n, e) is a normalization parameter chosen to match initial data against the 
Gaussian profile. Our result shows that for the k th order phase space Gaussian beam 
superposition, the following estimate holds on any bounded time interval, \t\ < T, 

(1.5) \K^-mr)\\iP<\\m^)-M\\^+\m\^- 

Here and in what follows we use A < B to denote the estimate A < CB for a constant C 
which is independent of e. 

For the initial data of the form ip[ n = A- m {x)e lSin<yX ^ e we need a superposition over an 
n-dimensional submanifold of phase space. The asymptotic solution is then represented 
as 

(1.6) i/> e (t, y) = Z(n, e) / iJ PGB (.t, y, X)6(w(t, X))dX, 

JQ(t) 

where w is obtained from the Liouville equation 

dtw + Hp ■ V x w — H x ■ VpW = 0, w(0, X) = p — V x Si n (x). 
Our result shows that 

(1-7) ||(^ e - ^)(i, -)|| i2 < |supp(A in )|et-f . 

We now conclude this section by outlining the rest of this paper: in Section 2 we start 
with Gaussian beam solutions in physical space, and define the phase space based GB 
ansatz through the Hamiltonian map. Section 3 is devoted to a recovery scheme through 
superpositions over a moving domain. The total error is shown bounded by an initial error 
and the evolution error of order e 1 / 2- "/ 4 . Control of initial error is discussed in Section 

4, followed by the convergence rate obtained for first order GB solutions. In Section 

5, we discuss how to use caustic structure to obtain some better error estimates. Both 
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convergence and convergence rate are obtained for higher order GB solutions in Section 
6. In Section 7, we present a new level set approach for construction of the phase and 
their derivatives. Finally in the appendix, we derive phase space equations for all involved 
Gaussian beam components. 

2. Phase space based Gaussian beam Ansatz 

2.1. First order Gaussian beam solutions. As is well known, the idea underlying 
Gaussian beams (GB) is to build asymptotic solutions concentrated on a single curve 
in physical space (t,x) € K x M n . This means that, given a curve 7 parameterized by 
x = x(t), one makes the Ansatz 

(2.1) i> e {t,x) =A € {t,x)e l ^ t ' x)/e , 

where <3?(t, x(t)) is real, and Im{&(t, x)} > for x ^ x(t). The amplitude is allowed to be 
complex and have an asymptotic expansion in terms of e: 

A e (t, x) = A (t, x) + eA^t, x) + • • • + e N A N (t, x). 

We wish to build asymptotic solutions to Pip(t,x) = 0, i.e., we want Pip e = 0(e 2 ). 
Substituting from (|2.ip . 

p^e = e iHt,x)/e ^ + H ( X} g x ^ Ao + e (_iLA + (dt$ + H{x, d x $))Ai)] + 0(e 2 ), 

where L is a linear differential operator, whose form is clear from f|2.3|) below. The key 
step in the GB construction is the choice of <I> such that + H (x, d x &) vanishes to high 
order on 7. The propagation of amplitude can then be determined by LAq = to make 
the 0(e) term vanish. We denote <3?(i,x) on the curve 7 by S and the leading amplitude 
by A = Aq. These lead to the standard WKB system 

(2.2) 8 t S + H(x, S/ X S) = 0, 

(2.3) d t A + H p ■ V X A = -^[Tr(H xp ) + Tr{V%SH w )], 

where Tr is the usual trace map. We then compute the Taylor series of dt<& + H(x, & x ) 
about x(t) to first and second order to obtain equations for the phase gradient and Hessian 
(u,M) = (V x $, V^$) as follows 

(2.4) d t u + H p -V x u = -H x , 

(2.5) d t M + H p ■ V X M + H xx + H xp M + MH px + MH pp M = 0. 

It is shown in [28] that the above construction is only possible if (x(t),p(t)), where p(t) = 
V x &(t,x(t)), is a (null) bi-characteristic curve, which is consistent with the characteristic 
system for the Hamilton- Jacobi equation <\2.2\) . 



(2.6) jx = H p , x(0) = so, 

(2.7) ± p =-H x , p(0)=po, 

where (x,p) = (x,p)(t; xq,pq). From here on we denote the phase space variable as X = 
(x,p) and Xq = (xo,po). Then the Hamiltonian dynamics can be expressed as 

(2.8) j t X(t, Xq) = V(X(t, X )), X(0, Xq) = X . 

The phase velocity V = (H p ,—H x ) is divergence free, i.e. divxiV) = 0. On this curve 
X = X(t, Xq), Gaussian beam components of the first order such as the phase S(t; Xq), the 



(2.9) -S(t;X ) = p-H p -H(x,p), 5(0; X ) = 5 in (x ), 



(2.11) -A(t; X ) = -- [Tr[H xp ] + Tr[MH pp }] , A(0; X ) = A in (x ). 
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Hessian M(t;Xo) as well as the amplitude A(t;Xo), are obtained by solving the following 
system of ODEs 

d 

~dt' 

(2.10) j t M(t; X ) + H xx + H xp M + MH px + MH pp M = 0, M(0; X ) = M in (X ), 

^(i;X ) = -- 
dt y 1 2 

The essential idea behind the GB method is to choose some complex Hessian M- m initially 
so that M remains bounded for all time, and its imaginary part is positive definite. This 
way the amplitude A(t;Xo) is ensured to be also globally bounded from solving (|2,lip . 
With these components in place, the Gaussian beam phase is constructed as 

$(t, y; X ) = S(t; X Q ) + p(t, X )(y - x(t, X )) + \{y - x{t, X )) T M(t; X )(y - x(t, X )), 

where p(t, Xq) = V x $(i, x(t, -X"o))- The leading order of the amplitude is taken as 

A(jt,y;Xo)=A(jt;X ). 

The above construction ensures that the following GB Ansatz is an approximate solution 

^ GB (t,y;Xo) = A(t;X )exp (l$(t,y;X ] 

The requirement that Im(M) be positive definite ensures that the asymptotic solution is 
concentrated on y = x(t,Xo), see e.g. [28] . 

2.2. Phase space based Gaussian beam ansatz. If we regard Xq to be the Lagrangian 
particle marker, then the map 

(2.12) X = X(t,X ) 

serves as a particle trajectory mapping: an initial domain Q £ R 2n in phase space evolves 
in time to 

(2.13) x(t 1 n) = {x(t,x ), i efi), 

with the vector V = (H p , —H x ) tangent to the particle trajectory in phase space. Since the 
velocity field is divergent-free, the elementary properties of X(t, Xo) tells that Vol(X(t, CI)) - 
Vol(n) = \fl\ and 

In other words the map is volume-preserving and invertible. 

The phase space based Gaussian beam ansatz is thus obtained by changing Xq to X 
through this particle-trajectory map: 

(2.14) i>PGB(t, y, X) = A(t, X) exp (U(t, y, X)J , 
where 

(2.15) y, X) = S(t, X)+p-(y-x) + ^(y- x) T M(t, X)(y - x). 

To derive the corresponding dynamics for (S, M, A) we need the following fact. 
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Lemma 2.1. ('Operator lifting ) Let the phase representative of w(t; Xq) be w(t,X) in 
the sense that w(t;Xo) = w(t, X(t, Xq) for any t > 0, then 

jw{t-X Q ) = Cw(t,X), 

where C is the usual Liouville operator defined by 

(2.16) C:=d t + V-V x . 
Proof. Taking differentiation of 

w{t; Xq) = w(t, X(t, Xq)), \ft>0 

in time, we obtain 

^-w(t; X Q ) = d t w + ^-X(t, X Q ) ■ V x w = d t w + V- V x w. 
dt dt 

□ 

Changing the time derivative ^ to the Liouville operator £ in the Lagrangian formula- 
tion of equations for (S, M, A) in (f279]) - ff27TT|) . we obtain the PDEs for (S, M, A) in phase 
space: 

(2.17) C(S)=p-H p -H(x,p), S{0,X)=S in (x), 

(2.18) C(M) + H xx + H xp M + MH px + MH pp M = 0, M(0, X) = M in (X), 

A(0,X)=A in (X). 

A method for solving the equation for M based on Ricatti equations was given by Leung 
and Qian (see [2], formulas (56)-(57)) . Jin, Wu and Yang have an alternative way of 
computing M based on complex level set functions (see [12J, formulas (3.5)-(3.16)). We 
give yet another method of constructing M, and the higher derivatives of the phase by 
level set methods in §7. We point out that though ipGB(t,y, Xq) is an asymptotic solution 
to the wave equation, ipPGB{t,y,X) is usually not. But we shall show that its integral 
over the moving domain X(t, r2(0)) remains an asymptotic solution of the wave equation. 
It is this remarkable feature that allows us to globally recover the original wave field from 
only some phase space based measurements! 



(2.19) 



C(A) 



2" 



Tr[H xp ]+Tr[MH pp ] 



3. Recovery of wave fields by superposition 

Since the wave equation we consider is linear, the high frequency wave field ip at (t, y) 
in physical space is expected to be generated by a superposition of neighboring Gaussian 
beams 



(3.1) r(t,y) = Z(n,e) ^ GB (t,y,X )dX , 

Jn(o) 

where 

Q(0) = {X , xq e supp(^i n ), p G range(d x S in (x))} 

is an open domain in phase space from which we construct initial Gaussian beams from 
the given data. The normalization parameter Z(n, e) is determined by matching the initial 
data ipo(y) so that 

l^oO-^o,-)!! -0, e^O. 
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^ e (i, y) = Z(n, e) / Vpgb(*, y, X)dX, 
Jn(t) 



By invoking the volume preserving map X = X(t,Xo) and its inverse Xq = Xq(1;,X), 
we obtain a phase space based Gaussian beam ansatz 

4>PGB{t,y,X) :=iJjGB(t,y,Xa(t,X)). 

As remarked earlier, since the map is time dependent, the above phase space based GB 
ansatz is no longer an asymptotic solution of the wave equation. We note that their 
superposition over the moving domain X(t, 12(0)) remains a correct asymptotic solution. 

(3.2) 
where 

n(t) = X(t,n(Q)). 

This can be seen directly by using change of variables to go back to the Lagrangian 
superposition (|3.ip . 

In what follows we will construct iPpgb without reference to i/jgb- While we could 
still recover ipGB from i/jpgb by coordinate transformations, we can check directly that 
superpositions of iPpgb are asymptotic solutions. This only requires the following two 
lemmas. 

Lemma 3.1. For any smooth f(t,X) and divergence-free velocity field V, one has 

(3.3) -| / f(t, X)dX = [ [d t f + V X ■ (fV)]dX. 

at Jx(t,Q) Jx(t,n) 

Our estimates are consequences of the following elementary lemma. 

Lemma 3.2. Assume that Im(&(t,y, X)) > c\y — x\ 2 , c > 0, and the Lebeque measure of 
the initial domain |^(0)| is bounded. Let B(t,y,X) be a smooth function, satisfying 



\B\ < C\y-x\ k , k>0. 



Then we have 



fi(t) 



B(t,y,X)e i ^ x)/e dX 



< 



n(0)|e 



k I n 
2 + 4, 



Proof. Using Minkowski's integral inequality we have 



fi(t) 



B(t,y,X)e i ^ t ' y ' X ^ t dX 



< 



\B\e~ Im ^l e dX 



n(t) 



1/2 



< 



n(t) \Jy 



\ B \2 e -2Im^)/e dy ^ 



dy 



1/2 



I// - x\ 2k e- 2c \y- x \ 2 l*dy)' 2 dX, 



< c 

ltl(t) \Jy 

continuing the estimate with the stretched coordinates y — x = e 1//2 y', and changing from 
y to y' in the integral 

1/2 



< C [ £ f + f ( [ \y'\ 2k e- 2c \y'\ 2 dy>) dX 
Jn(t) \Jy' J 



C\n(t)\e 



1 T 4 



2k *- 2c \y\ 2 d y 



\y\ e 



1/2 
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which when using |0(t)| = |^(0)| proves the result. □ 

The normalization parameter needs to be chosen to match the initial data. For example, 
if initially Im(M- ia ) = (51, (5 > 0, then we need to arrange to match the initial data against 
exp(— (3\x — y\ 2 /e). That accounts for 

y T Im(M in )y/(2e), \ 1 _ ( _/3_V /2 n /2 



(3.4) Z(n,e)= (J e -V '^JWW^j = ^ e 

in dimension n. Taking the Schrodinger equation as an example, we obtain the following. 

Theorem 3.3. Let P be the linear Schrodinger wave operator of the form P = —iedt + 
H(y, —iedy), where H(y,p) = — h V ex t(y), and ip e is defined in h3. 2i) with Im(M) being 
positive definite, and Z(n,e) ~ e~ n / 2 ; then tfj e is an asymptotic solution and satisfies 

(3-5) \\PW)(t,-)\\ L ><\n(0)\e*~. 
Proof. We apply the operator P to both sides of (|3.2|) to obtain 

Z- X PW\ = (-ied t + H(y, -ied y )) [ $ PGB (t, y, X)dX 

JQ(t) 

(3.6) = l [P[iPp GB }-ieVx-(Vip PGB ))}dX. 

JU(t) 

By a straightforward calculation it follows 

P^pgb] = -ied t A(t,X)e^ (t ' y ' x)/e + AP[^ y ^l e ] 



A[d t $ + H(y, d y $)] - ie ( 8 t A + ^Tr[dfe] 



The transport term in the integrand gives 

-teV x ■ (V^ PGB ) = e *('**)A 



-ieV ■ V X A + AV ■ V x $ 



Putting together we have 
(3.7) 

A, 



[P-ieV x -V](*pPGB) = e 



i$(t,y,X)/e 



A[C[$]+H(y,d y $)] - ie C[A] + -Tr[^$] 



Using = M(t,X) and that A solves (I2J9D . we see that 0(e) term vanishes. From 
(j2TT5j) it follows that 



and 



C[^\ = C[S] - \p\ 2 - d x V ext ■{y-x)+ l -{y - x) T C[M](y - x)-p T M(y - x) 
H(y,dy$) = V ext (y) + ^\p + M(y-x)\ 2 . 



This when using equation (|2.17p for S gives 
(3.8) 



C[$] + H(y,d y $) = V ext (y)-V ext (x)-d x V ext (y-x) + ^(y-x) T d 2 x V ext (y-x) = 0(\y-x\ 3 ). 



1 

+ 

Consequently, using the above lemma with k = 3, 

||P[^](t,-)llL ? <^, e )l^(0)|6l + ?, 

which with (13.41) leads to the desired estimate. □ 
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Remark 3.1. The extra term in the integral in (|3,6p gives an alternate way of seeing that 
the phase space super-position is an accurate solution of the PDE. Of course, it integrates 
to zero when the support of the beam superposition does not touch the boundary of the 
integration domain, but it makes it possible to verify the accuracy without going back to 
the Lagrangian super-position. 

We now obtain the following estimate. 

Theorem 3.4. Given T > 0, and let ip be the solution of the Schrddinger equation subject 
to the initial data ipo, and vp € be the approximation defined in \3. ty) with Im{M^ n ) being 
positive definite, and |^(0)| < oo. Then there exists eo > 0, a normalization parameter 
Z(n, e) ~ e~ n / 2 , and a constant C such that for all e G (0, eo) 



- t/,)(t, OIU* < 11^(0, •) - + C|O(0)|e2-4 

forte[0,T]. 

Proof. Let e := ip e — if), then from P[ip] =0 it follows 

P[e] = P[^ e } - P[ip] = P[^ e ]. 
A calculation of J R „[eP[e] — eP[e]]dy leads to 

e— J \e\ 2 dy = J Im(eP[Tp e ])dy. 



Integration over [0, t] gives 



1 



(3.9) Mt,-)\\iZ<M0,.)\\ L *+-J \\Pm(r r )\\ L 2dr, te[0,T\. 

This when combined with the estimate for -PfV' 6 ] in (|3.5p gives the result as desired. □ 

Remark 3.2. The approximation error comes from two sources: initial error and the evolu- 
tion error. To improve accuracy one has to enhance the accuracy for both. The evolution 
accuracy can be improved by obtaining more phase space based measurements such as 
higher order derivatives of phase and amplitude, which will be sketched in §5. 

Remark 3.3. In phase space the tracking of the beam propagation is lost in the Gaussian 
beam ansatz, but has been recorded through the moving domain f2(i), which can be traced 
back to n(0). 

4. Control of the initial error 

Let K(x,t) = ^ nT ^ n / 2 e At be the usual heat kernel, satisfying \\mit T ^K{x, r) = 6(x) 
as distributions on W 1 . Then 

I K(x - y, r)dx = 1, Vr > 0, y G M n . 

For highly oscillatory initial data we have 

^in(y) = AUy)e iS ^ y)/t = j A in (y)e iSi ^K (x - y, |) dx. 

Both the phase and amplitude in the integrand can be approximated by their Taylor expan- 
sion when \x — y\ is small, say \x — y\ < e 1 / 3 , and the integral will then be 0(exp(— c/e 1 / 3 )) 



10 HAILIANG LIU AND JAMES RALSTON 

with some c < | outside this neighborhood. Let Tj[f](y) denote the j th order Taylor 
polynomial of / about x at the point y. Then 

(4.1) Vin(y) ~ J AUx)eil T ^ Si ^K (x - y, |) dx, 

which tends to ip\ n as e — > 0. 

Indeed, the approximate accuracy is ensured by the following result by Tanushev [33J. 

Lemma 4.1. Let S- m £ C°°(R n ) be a real-valued function, and A- m 6 CQ°(R n ), and 
p £ C^°(R n ) be such that p > 0, p = 1 in a ball of radius 5 > about the origin. Define, 

^n{y) = AUy)e lSUv)l \ 

e 



;(y;x) = p[y - x)Tf[A, n ]{y)e^ T ^ s ^K (x-y, 2 



Then 



V'in(-) - / v(-;x)dx 

supp(A in ) 



L 2 



< e 2 , 



Remark 4.1. We note that for the case j = the cutoff function p is unnecessary, and 
it can be shown the initial approximation error remains of order e 1 / 2 . However, a cutoff 
function is certainly important when one is building beams of higher accuracy because 
the higher order terms in the Taylor expansion of the phase can change the sign of its 
imaginary part when one does not stay close to the central ray, see §5. 

We take the above approximation in (14. lj) as initial data for ip e (0,y) and rewrite it as 
follows 

(4.2) ^(0,y) = Z(n,e) [ ^ PGB (0,y,X)5(p - V x S in (x))dX, Z(n,e) = \ 

with Ai n = Ai n (x), Si n (X) = Sin(x) and M\ n (X) = d^.S\ n {x) + il. The above lemma 
ensures that 

(4-3) ll^in(-)-^(0,-)llL2<e 1/2 . 

In order to track the deformation of the surface p — V x .5 , ; n (x) = as time evolves, we 
introduce a function w = w(t, X) such that 

(4.4) C[w) = 0, w(0, X)=p- V x S in (x). 

For smooth Hamiltonian H(x,p), w remains smooth once it is initially so. A modified 
approximation is defined as 



(4.5) V £ (i, y) ■= Z(n, e) / ^ PG B(t, y, X)5(w(t, X))dX, 

JQ(t) 

which has taken care of the Dirac delta function in (|4.2p . We then have the following 
theorem. 



Theorem 4.2. // assumptions of Theorem 13.31 are met, then ^ defined in fl^.5[ ) is also 
an asymptotic solution, satisfying 

3 71 

(4-6) \\P[r](t,-)]\\ L 2 < |supp(A m )|e2-4. 
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Proof. Using the volume-preserving map of X = X(t,Xo) and w(t, X(t, Xq) = w(0,Xo), 
we have 



^(t,y) = Z(n,e) / i/> PGB (t,y,X(t,X ))5(w(t,X(t,X )))dX 
IQ(0) 



= Z(n,e) / ^ GB (t,y;X )5(w(0,X ))dX 
Ja(o) 

= Z(n,e) / xfj G B(t,y;X )d(po - V x S in (x ))dX 
(4.7) =Z(n,e)[ A(t;x )G(t,y;xo)e^ x ^dx . 

Jsupp(A in ) 

Here for simplicity we use only xo in the integrand instead of Xq = (xo,V x Si n (xo)). We 
now repeat the similar estimate to that in the proof of Lemma 13.21 with k = 3 to obtain 



(4.8) 
Hence 



A(t;x )G(t,y-x )e l ^y^^dx 

supp(A in ) 



< !supp(^ in )|e2 + 2. 



3 i n 



\\P[^](t,-)\\ L 2 <Z(n, e )j|supp(A in )|e2 + 4 

which with (|3.4p leads to the desired estimate. □ 

Plugging estimates (|4.3p and (|4.6p into (|3.9p we arrive at our main result. 

Theorem 4.3. Given T > 0, and let ip be the solution of the Schrddinger equation subject 
to the initial data ifj- m = A- m e lSin ^ x ^ e , and ip € be the first order approximation defined in 
J^.5[ j with initial data satisfying Si n (X) = Si n (x), M- m {X) = d x S m (x) +il, and A- m (X) = 
Ai n (x) with \supp(A- m )\ < oo. Then there exists eo > 0, a normalization parameter Z(n, e), 
and a constant C such that for all e € (0, eo) 



-VOfoOlU 3 ^ |supp(A in )|e2 4 
forte [0,T]. 

Remark 4.2. The exponent 1/2 reflects the accuracy of the Gaussian beam in solving the 
PDE. It will increase when one uses more accurate beams. The exponent — f indicates 
the blow-up rate for the worst possible case due to caustics. Of course, if nature of the 
caustic were a priori known, it would be possible to obtain a better convergence rate by 
taking the caustic structure into account. 

5. A CLOSER LOOK AT CAUSTICS 

5.1. Schur's lemma. Instead of using the Minkowski inequality we shall use Schur's 
lemma to see how caustic structure may be used to obtain a better estimate. Recall 
Schur's Lemma: If [Tf](y) = J K(x,y)f(x)dx and 

su Pz / \K(x,y)\dy = d, sup y / \K(x, y)\dx = C 2 , 

J y J x 

then 



|T/|| L2 < s/C&\\f\\v. 
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Proof. We have by Schwartz 

\[Tf](y)\ 2 < (I \K(x,y)\f(x)dx) < I \K(x,y)\dx I \K(x,y)\\f(x)\' 2 dx 



<C 2 J \K(x,y)\\f(x)\ z dx. 

So integrating both sides in y and taking the square root gives the result. □ 
We now apply Schur's lemma to left hand side of (|4.8j) . So for simplicity 



[Tf](y)= / A(t;x )G(t,y;x )e i ^ x ^dx , 

Jsupp(A in ) 

where the imaginary part of $>(t,y;xo) is bounded below by cl and for convenience we 
will assume that \G\ < \y - x(t,x )\ k . Then one can apply Schur's lemma with 

d = sup X0 [ \y - x(t,x )\ k e- {c/e)ly - xit ' Xo) f 'dy = et+f f \z\ k e~ c ^ 'dz, and 



C 2 (t,e) = su Py / \y-x(t,x )\ k e-^- x ^\ dx . 

In general one does not know what C^if, e) will be. As long as A has compact support 
C 2 will be at least bounded by ce fc//2 . Thus the error in L? norm will be bounded by 
ce fc//2+n / 4 , as shown in (|4.8p with the Minkowski inequality. Note that the worst case for 
the wave equation, Ofy — Aip = 0, occurs when x(t, xq) = xq{1 — t/\xQ\). In that case 
C2 = ce fc / 2+ ^ 1+n ) //4 , yielding a better rate. 

5.2. An example with remarkable accuracy. Here below we illustrate that a bet- 
ter convergence rate can also be obtained for the Schrddinger equation with quadratic 
potential and quadratic phase. We consider the solution of 

ied t t/) = -—Aip, x G M n , 
with the initial data ip(0,x) = exp(— i|x| 2 /(2e)), then 

^(i,*) = (i-r«/ 2 ex P (- 2e(1 _ t) 

This solution becomes a multiple of the 5- function at t = 1. This suggested solving the free 
Schrddinger equation with initial data ip(x,0) = g(x)exjp (— i|x| 2 /(2e)) where g S C°°(R n ) 
and evaluating the solution at t = 1. Using Fourier transform one may express the solution 

Evaluating at t = 1, we have 

= {2 J e)n /2 j y 9{y)e- lx - y/e dye^\ 2 /^ = ce^gix/e)^^^ , 

where c = e -7rm / 4 and g is the Fourier transform of g, defined by 

1 

(2tt)™/ 2 

It is easy to verify that 



9(0 = tA^ I g{x)e**dx. 
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As e — > 0, ip(l, x) diverges (pointwise) like e~ n l 2 near x = 0, but goes rapidly to zero away 
from x = 0. 

We now build a superposition of Gaussian beams approximation for the solution of the 
same problem. For the Gaussian beam superposition, the phase is obtained as 



Ht,x;y) = (t-l)\f-x-y + \y\M-t) , 1 + {fJ ._ 1)t 2 • 

Note that we have chosen (3 > for the initial beam width. For the amplitude we get 

(1 + (f3i - l)t) n / 2 A(t, y(l - t)) = A(0, y) = g(y). 
So, setting x(t;y) = (1 — t)y, we obtain 

A{t,x{t-y)) = (l + (f3i-l)t)- n l 2 g(y). 
If we do the superposition with the normalization, we end up with 



^(t,x) 



V2vre 



[l + (J3i-l)t]- n ' 2 g(y)e 



J$(t,x;y)/e 



dy. 



If we evaluate that at t = 1, it becomes 

n/2 



V27re 



[(3i\- n ' 2 g{y)e^- x - y+§ Tr^dy 



= ce- n ' 2 g(x/e)e^ 1 ^ 2 /^ 
= i>(l,x)e~W 2 ^. 

This shows that at the caustic x = 0, both become the same. We can see the error 
ip(l,x) (l — e - ^! 2 '! '( 2/3e )^ when measured in L 2 -norm: 

||^(1, •) - m OHi, = 6- / \g ) | 2 (l - e -NV(W) 2 dx 
= f \g(z)\ 2 {l-e-^l^) 2 dz. 



That implies that 

(a) For any g G 1? the Gaussian beam approximation converges to the true solution (at 
t = 1), but there is no uniform estimate on the difference in terms of the L 2 -norm of g (an 
example of strong but not uniform convergence). 

(b) If J \g(z)\ 2 (l + \z\ 2 ) 2 dz < oo, i.e. if g £ H 2 , then the norm of the difference is 0(e). 
Actually in the current example, it can be verified that the evolution error is zero, so 

the initial error should propagate in time. If we look at the initial error of the Gaussian 
beam approximation, we have 



(0,-)-^(0,-)|li a =j[ (j^) f y9 (y)e-^y\ 2 dy-g{x) 
= i^l) 1 f x \f(9(y)- 9 (x))e~^ x ^ 2 dy 



dx 



dx. 
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Set K{x) = (^J e'^ lxl , a direct integration shows that K = (2TT)- n / 2 e~ e ^ 2 ^ 2 ^ . If 
we apply Parseval's theorem, we obtain 

11^(0, ■) - m-)\\h = B - (^) n/2 g ■ Kf L2 

(5.1) = / \Hz)\ 2 (l-e-M 2 'W>) 2 dz, 



which is the same error as that evaluated at t = 1. We point out that for the initial phase 
of general form, the initial error is still 0(e 1//2 ) unless an higher order expansion of the 
phase is used. 

There are two conclusions that one can draw from the preceding. 



Theorem 5.1. Under assumptions of Theorem \4-3\ end assume that the potential is a 
quadratic function. Then for t G [0, T] and e G (0, eo) we have 



• If Sin ^ a quadratic function and A\ n £ H 2 



• If Si n e C°° and An e C °° 

\W-mr)\\^<e l/2 . 
Proof. It follows from (|3.8p that for quadratic potentials 

P[ip e ] = 0. 

Then the total error is governed by the initial error only. For quadratic potentials and 
A in € H 2 (R n ) we obtain the O(e) error as shown in (|5,ip . For the general phase function 
the claim follows from Lemma 14. II with j = 0. □ 

6. Higher order Approximations 

The accuracy of the phase space based Gaussian beam superposition also depends on 
accuracy of the individual Gaussian beam ansatz. If we refer the above construction as a 
first order GB solution, then the k th order GB solution should involve terms up to (k + l) th 
order for the phase, and (k — 1 — 2l) th order for the I th amplitude Ai for I = 0, • • • , [— ip] ■ 
The equations for these phase and amplitude Taylor coefficients can be derived by letting 
the leading order ones to hold on 7 along with several of their derivatives (see the Appendix 
)• 

Let X = X(t; Xq), with x = x(t; Xq), denote the characteristic path at time t > 0, which 
originates from Xq. Following [33] we define the k th order Gaussian beams as follows 



ipkGB(t,y;X ) = p(y - x) 



1=0 



exp(^+i[*](y) 



where T?[f](y) is the k th order Taylor polynomial of / about x evaluated at y, and p is 
a cut-off function such that on its support the Taylor expansion of $ still has a positive 
imaginary part. 

By invoking the volume preserving map X = X(t,Xo) and its inverse map denoted by 
Xq = Xo(t,X), we obtain a phase space based k th order Gaussian beam Ansatz 

^kPGB(t, y, X) := il) kGB {t, y, X (t, X)). 

Beyond the first order GB components, all Taylor coefficients dy $ for \a\ > 3 in phase 
space are replaced by m a (t, X), satisfying a linear equation (|8.ip in phase space; and 
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Taylor coefficients dyAi for \a\ > 1 in the amplitude are replaced by Ai j0l , which can be 
obtained recursively from solving transport equations in phase space, see the appendix for 
details. 

Proceeding as previously, we form the superpositions. 

(6.1) ip e k (t, y) = Z(n, e) I ^ kPGB {t, y, X)6(w(t, X))dX, 

J(i(t) 

where Q(t) = X(t,£l(0)), and w(t,X) is the solution of the Liouville equation subject to 
w(0,X)=p-V x S in (x). 

This gives a k order asymptotic solution of the wave equation. More precisely, we 
have the following theorem. 

Theorem 6.1. Let P be the linear Schrodinger wave operator of the form P = —iedt + 
H(y, —iedy), where H(y,p) = ^ — h V ex t(y), and ip e is defined in 116. 1\) with Im(M Ui ) = I 
and Z(n,e) = (2ne)~ n l 2 , then tf)f is an asymptotic solution and satisfies 

(6-2) ||P[^](t,-)ll^<|supp(A in )| e l +1 -?. 

Proof. Using the volume-preserving map of X = X(t,Xo) and w(t, X(t, Xq)) = w(0,Xq), 
we have 



V>|(i, y) = Z(n, e) / ^ kPG B(t, V, X(t, X ))5(w(t, X(t, X )))dX 
'n(o) 



Z(n,e) / ^ kGB (t,y;X )6{w{0,X ))dX 
1(1(0) 



Z(n, e) / ip kGB (t, y; X Q )5(p - V x S in (x ))dX 
J(i(o) 



= Z(n,e) ifj kGB (t,y;x )dx . 

J supp{A in } 

According to the GB construction sketched in the appendix, tpkGB(t, y\ xo) are asymptotic 
solutions for each xq, so will be their superpositions il) e k {t,y). It remains to verify (|6.2|) . 
First we see that 

P[^%{t,y)] = Z{n,e) I P[ip kGB (t,y;x )]dx . 



Z(n,e) [ P 

Jsupp{A in } 



EU 1 z l TU-2i[M{y) 



and $ 



Using (|8.3p in the appendix with A replaced by p(y — x) 
by T% +1 [$](y), we have 

Co (t,y) = [d t T£ + M(y) + H(y,v y TZ + Mmp(y-x)TLMo](y)- 

Using T% +1 [$>](y) = 5>(y) + Rf +l [$>](y), here R k+1 denotes the remainder of the Taylor 
expansion, and G(t, y) = dt$ + H(y, V y $) = 0(\y — x\ k+2 ) we can see that 



\co(t,y)\ <C\y - x 



fc+2 



Also using the construction for A\ and their derivatives in the appendix, we are able to 
show 

\ Cl (t,y)\<C\y-x\ k+2 - 21 , 
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where we have used the fact that differentiation of p vanishes in a neighborhood of x. 
The use of the cut-off function ensures that we can always choose a small neighborhood 
of x(t; xq) so that 

Jm(T£ +1 [$](y)) >c\y-x\ 2 . 
Consequently, using Minkowski's integral inequality, 

Z- 1 \\Pmt,-)]\\L*< \ I I e -^ +1 [*Ky))A| C0 + Cie + ...| dxc ' 

'supp(An) 





supp(A in ) \Jy 




1/2 

e -2c\y-x(t,x )\ye ^ + ^ + . . . |2 dy \ ^ 



<c/ ( [ e -Mv-*<t,*o)\V* |y-*(i,x )| 2(fc+2 ~ 2 V^ dx . 

Jsupp(A in ) \Jy l=0 I 

If we introduce the stretched coordinates y — x(t; xq) = e l / 2 y' , and changing from y to y' 
in the integral, we see that the new integrand is bounded by 

e fc+2+f| y ,|2(fc+2-20 exp (-2c|y'| 2 ). 

Thus •)] IIl 2 is bounded by Z(n, e)|supp(Ai n )|e2 +1+ 4 . The desired estimate then 

follows. □ 

In order to obtain an estimate of \\{ip% — tl?)(t,-)\\ for any t < T, all that remains to 
verify is that the superposition (|6.ip accurately approximates the initial data. For t = 0, 
the approximation is as follows 



#(0, y) = Z(n, e) / Apgb(0, y, X)5(w(0, X))dX 
Jn(o) 



Z(n,e) / ip kG B(0,y,X )5(po -V x S in (x Q )))dX 

JQ(0) 



= Z(n,e) / ip kGB (0,y,xo,V x Si n (xo))dxo, 

J SUppAi n 

where 

i> kGB (0,y,x,V x S in (x)) = p(y - x) exp (~2f +1 [Sy (y)) e'^ 2 ^, 

where we have taken A = A in and Ai = for I > 1, d%$(0,x) = d°S- m (x)(a ^ 2), and 
d 2 3>(0, x) = d 2 S- m (x) + il. From Lemma 14. II we have that 

Uin-^(0,-)\\ L 2<el 
Thus our main result for k th order phase space GB superposition is as follows. 

Theorem 6.2. Given T > 0, and let ip be the solution of the Schrddinger equation subject 
to the initial data ip{ n = A{ n e iSin ^ x ^ e , and ^ e be the k th order approximation defined in 
with initial data chosen as described above with \supp(Ai n )\ < oo. Then there exists 
eo > 0, a normalization parameter Z(n,e) ~ e -n / 2 , and a constant C such that for all 
e£(0,e ) 

\W-m,-)\\L* < |supp(A in )|ef-t 

for t € [0, T]. 
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7. Computing Taylor coefficients of the phase via level set functions 

We now turn to construction of the phase space ingredients required for the approx- 
imation. In order to identify a bi-characteristic curve in phase space, we introduce a 
vector- valued level set function <f> E M. 2n so that the interaction of zeros of each component 
uniquely defines the target curve. In other words, we assume that 

T = {(t,X), <j ) (t,X) = ( f>(0,X )} 

contains the bi-characteristic curve starting from Xq = (xo>Po) f° r anv t > 0, then <p must 
satisfy 

4>(t,X(t,X )) = </>{0,X Q ). 
This is equivalent to the following Liouville equation 

(7.1) C[<Kt,X)]=0, 

where C := dt + V ■ Vx is the Liouville operator. The initial data can be simply taken as 

(7.2) </>(0,X)=X-X o . 

Then the curve T is globally determined by the zero set of a vector level set function 
= (0i,0 2 ) T . 

For the construction, S can be solved from (|2.17p . we are then left to determine M, 
followed by solving (|2.19|) to obtain A. Note that equation (|2.18|) is nonlinear in M, the 
solution might not exist for all t > 0. The heart of the GB method is to choose complex 
initial data so that a global solution is guaranteed and satisfies two requirements [28J : 
i) M = M T , 

\\)Im{M) must be positive definite for all t > 0. 

7.1. Evaluation of the Hessian. We now show this can be done via the obtained level 
set functions 4> € M. 2n . 

Theorem 7.1. Let (ft = (^i,^2) T with G R n be the global solution of \7. 1\ ) with the 
initial condition |7.^| ). We have 

a) C(k\(f>i + k2(p2) = for any k\, k-2 € C. 

b) Set g := k\4>\ + ki4>i- If Imykife) ^ 0, then g p is invertible for all t > 0. 

c) If M- m = -gxigp)' 1 ^, then M = -g x {g p )~ l for all t > 0, 

d) M satisfies 12.18\) and i). If Im(k\/k2) < 0, then M satisfies ii) too. 

Proof, a) This follows by noting that the Liouville operator is linear and all its coefficients 
are real. 

b) By taking the gradients V x and V p of the Liouville equation C(g) = 0, respectively, we 
obtain the following equations 

(7.3) £{9x) = Hxxdp — Hxp9x, 

(7.4) C{g p ) = H px g p - Hppg x . 

The equation is understood to be satisfied by each matrix. Let B = ~g p sr g x — g~ x r g p be a 
complex matrix, and / an identity matrix, a direct verification shows that 

C{B) = C{gfg x ) - C{gf g p ) 



£{g P ) 9x + 9 P C(g x ) - C{g x ) g p + g x C{g p 
0. 
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Observe that B(0,X) = —2iIm{kxk<z)I is a constant matrix. Thus for any t > 0, 
B(t,X(t,X )) =B(0,X Q ) = -2iIm(kik 2 )I- 

The condition Im(k\k2) 7^ ensures that g p must be invertible for all t > 0. Otherwise 
there would be a nonzero vector c such that g p c = 0, hence c T Be = (g p c) T g x c = 0, leading 
to a contradiction. 

c) Set Q = g x + Mg p . A calculation using (I2.18j) . (f773T) and fTil) gives 

C[Q) = C(g x ) + C(Mg p ) = -(H xp + MH pp )Q. 

If M in = -g x {g p )~ l initially, then Q(0,X) = for all X € M 2n . Thus we have 

Q(t,X) =g x + Mg p = 0. 

This gives 

M = -g x (g P )~ l 

for all t > since <7 P is invertible. 

d. i) Initially M in = -jj^I = M^. Since M T also satisfies equation (f2T8]) . hence M = M T . 
ii) With the definition of B, we have 

B = -g^ T Mg p + gfMg p = -2ilm\gf Mg p ]. 

Initially we have g v = k2l- This together with B(t, X) = B(0, Xq) along T gives Im[g p r Mg p ] = 
\k2\ 2 Im[Mi n ]. Note that Im[Mi n ] is positive definite, hence Im[M] remains positive defi- 
nite for all t > 0. □ 

Remark 7.1. The formula M = —epx^epp)^ 1 , first derived in [10], plays an important role 
in [10] in deriving the equation 

£[/] = o 

for the quantity f(t,X) = \A(t,X)\ 2 det((j) p ), which remains globally bounded even when 
(j) p becomes singular. We note that a complex level set function was used in [12] to obtain 
a globally bounded Hessian. 

Remark 7.2. Since the Liouville equation is geometric and homogeneous, for each fixed Xq, 
the shift Xq in the level set function can be simply ignored, and be added back whenever 
it is needed. In other words, we can take initial data (p(0, X) = X, then the curve T can 
be represented as Xq level set: 

T = {X, ( j ) (t,X) = X }. 

Remark 7.3. If we follow this construction, the initial data for Mi n then depends on 
how we initialize the level set function <j). If we take k\ = [3 > and k 2 = i, then 
g = (3cpi + i(p2- If ((j>i, (p2)(x, X) = (x,p), then M- m = i(3I. If po is restricted to be 
the phase gradient at xq initially, then the initial level set function can be chosen as 
(01 , <p2){0,X) = (x,p — y x S- m (x)), this leads to M- m = d x S\ n (x) + i/3I. In this case 4> 2 is 
the function w(t,X) from (|4.4p . 
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7.2. Evaluation of higher order derivatives of the phase. Let G(t, y) = dt& + 
H(y,Vy&). If one wants to have G(t,y) vanish to a higher order than two on 7, it is 
necessary to obtain higher order derivatives of <£. In the appendix we derive a system of 
linear equations for m a (t,X) = dy&(t, x(t, Xq)) on 7. We now show that this again can 
be done through the vector-valued level set function (p. 

Differentiating (pi(t, x, V x $) = 0, I = 1, 2, to order of r > 3 we obtain 



i=l \P\=r 

for all multi- indices a of length r. Let g = k\(f)i + /c2</>2, again using the invertibility of g p 
we can obtain 



V x (m c 



si L <*v 

\v\ =r 



We do this recursively, since the coefficients c Qjr) = (ci ai?? , ■ 



T and d r 



(diai ' ' ' d nct ) 



depend on all the partials up to order r — 1. Since g p is invertible, the obtained derivatives 
remain bounded for all t > 0. 



8. Appendix 

In this appendix we follow [29J to determine higher order derivatives of phase and 
amplitude on 7, and further derive phase space equations they satisfy. From dyG = on 
7 with |q| > 3, we obtain 

d t (d^) + H p -v x (d^)+ Y, c ajV d^ + d a = o, 

\r)\=\a\ 

where c a ^ and d a depends on dy $ for \k\ < \a\. Using the Hamiltonian equations = H p 
we obtain 

j t (d^(t,x(t,X );X )+ Y c a , v d^(t,x(t,X );X ) + d a = 

\v\=H 

on (t,x(t,Xo)). Following Lemma 12.11 we obtain a linear system of Liouville type PDEs 
for partial derivatives m a (t, X) of a fixed order: 

(8.1) C[m a ] + Y Ca^v + d a = 0. 

\v\=H 

We solve the system (|8.ip starting from \a\ = 3, then |a| = 4 and so on until \ot\ = k + 1 
for k th order GB solutions. Since equations are linear, we have solutions defined for all 
t > 0. This construction ensures that 

(8.2) G(t,y)=0(\y-x\ k + 2 ). 

To determine the Taylor series of Ai, I = 0, • • • , N on 7, one proceeds as follows. Define 
the coefficients Cj (t, y) by 

(N+2 \ N 

J>(t,y)eM e^, A = ^V- 
3=0 I j=0 
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Then, with P = —iedt + H(y, —ied y ), we obtain 

co(t,y) = G(t,y)A , 
ci(t,y) = -iLAq + G(t, y)A 1 
cj+i(t, y) = -iLA x + G(t, y)A l+1 + g h I = 1, • • • AT + 1, 

where L is a linear differential operator with coefficients depending on 

1 



d t + H p ■ V y + - [tr(#«p) + tr(MH pp (y, $„))] , 



and 5/ = --jAj^-i. 

Thus to make P[ip e ] = 0(e K ) for a given K € Z, we now only need to make Cj vanish on 
7 to sufficiently high order. To do so we can solve the equations LA\ + igi = recursively 
starting with 1 = (go = 0), and solve it to arbitrarily high order by solving the linear 
transport equations for the partial derivatives of A\ that one gets by differentiating the 
above equations. From the above procedure we see that the number of terms, N, in the 
solution ansatz for border GB approximation is determined by the following relation 

fc-1 jfe+1 
< N < 



In other words, N = L^pJ + 1- Actually, given G(t,y) vanishes to order k + 1 on 7, we 
can choose the Taylor series of Aq on 7 up to order k — 1 so that C\ vanishes to order k — 1 
on 7. Passing to the higher order equations q+i = 0, we see that we can choose A\ so 
that q + i vanishes on 7 to order k + 1 — 2(1 + 1). Thus we need 2 > k + 1 — 2A^ > 0. 



Thus for k th order GB solutions, it is necessary to compute d% A\ for \a\ < k 



21: 



L(d«A )+ J2 
hl<M 



a 
V 

a 
V 



d^L^Ao) 



= 0, \a\<k-l, 



d^L(d^A t ) - -Ayd^A^ 



0. 



\a\ < k- 1 -2L 1 = 1, 



Lifting the operator into the phase space, we can obtain A[(t,X) recursively by solving 

A, 



tr(H yp ) + tr(M(t,X)H pp (y,p)) 



9l, 



where C is the Liouville operator. Same lifting can be applied to all involved derivatives 
of the amplitude Ai for I = 0, • • • [^f^J • This completes the construction for all involved 
Taylor coefficients of both phase and amplitude. 
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